home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / r_check3du.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  5.4 KB  |  243 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3.  
  4. #include "fft.h"
  5. #include "constant.h"
  6.  
  7. /*                        */
  8. /*    Precision dependant declarations    */
  9. /*                        */
  10. #ifdef    DOUBLE
  11.  
  12. #define TOLERANCE   DTOLERANCE
  13. typedef        double    this_half;
  14. typedef        double    this_type;
  15.  
  16. #define        THIS_SUM    dsum_
  17. #define        THIS_GEN    dgen_
  18. #define        THIS_FTU    dft3du_
  19.  
  20. #define        THIS_FFTUI    dfft3dui
  21. #define        THIS_FFTU    dfft3du
  22. #define        THIS_SCAL    dscal3d
  23.  
  24. #endif
  25. #ifdef    SINGLE
  26.  
  27. typedef        float        this_half;
  28. typedef        float        this_type;
  29.  
  30. #define TOLERANCE   STOLERANCE
  31.  
  32. #define        THIS_SUM    ssum_
  33. #define        THIS_GEN    sgen_
  34. #define        THIS_FTU    sft3du_
  35.  
  36. #define        THIS_FFTUI    sfft3dui
  37. #define        THIS_FFTU    sfft3du
  38. #define        THIS_SCAL    sscal3d
  39.  
  40. #endif
  41.  
  42. /*                        */
  43.  
  44. this_half THIS_SUM();
  45.  
  46. void inimat_();
  47. void GetArguments();
  48. void get_values();
  49.  
  50. int size_x, ld1, ld2, size_z , size_y, n_trials;
  51. this_type *pa, *pb, *pRef, *pSave;
  52.  
  53. main(argc,argv)
  54. int argc;
  55. char *argv[];
  56. {
  57.     int i, j, k, n_total, is_wrong, nerr, inc;
  58.     double x, y, dx, dy, emax, sum, err, maxerr;
  59.     this_half    t;
  60.     this_type *ptr, *pt;
  61.  
  62. /* ******************************************************* */
  63.     GetArguments( argc, argv);
  64. /* ******************************************************* */
  65.  
  66.     srandom( (123*getpid()) | 0x01);
  67.  
  68.     for( ; n_trials > 0; n_trials --) { 
  69.     get_values();
  70.  
  71.     printf("\n");
  72.     printf( "%3d <%3d,%3d,%3d>", n_trials, size_x, size_y, size_z);
  73.     fflush(stdout);
  74.  
  75.     ld1 = ((size_x + 2)/2)*2;
  76.     ld2 = size_y + 1;
  77.     n_total = (ld1*ld2*(size_z+1));
  78.     pa = (this_type *)malloc
  79.         ( (3 * n_total + size_x + 3* size_y + 3*size_z + 10 * FACTOR_SPACE) * sizeof(this_type));
  80.     if( !(pa) ) {
  81.         fprintf( stderr, "Could not allocate ... Exiting\n");
  82.         exit (-1);
  83.     }
  84.     pb = pa + n_total;
  85.     pRef = pb + n_total;
  86.     pSave = pRef + n_total;
  87.  
  88. /* *******************************************************
  89.     Initialize arrays
  90. ******************************************************* */
  91.     THIS_GEN(pRef, &n_total);
  92.     bcopy( pRef, pa, n_total*sizeof(this_type));
  93.     bcopy( pRef, pb, n_total*sizeof(this_type));
  94.  
  95. /* *******************************************************
  96.     NOT SO PACKED direct Fourier Transform
  97. ******************************************************* */
  98.         j = -1;
  99.     THIS_FTU( &j, &size_x, &size_y, &size_z, pb, &ld1, &ld2);
  100.     pSave = THIS_FFTUI( size_x, size_y, size_z, pSave);
  101.     is_wrong = THIS_FFTU( -1, size_x, size_y, size_z, pa, ld1, ld2, pSave);
  102.  
  103.     emax = TOLERANCE*(size_x * size_y * size_z );
  104.     sum = 0.;
  105.     err= 0.;
  106.     nerr= 0;
  107.     maxerr= 0.;
  108.     for( k = 0; k < size_z ; k++)
  109.         for( j = 0 ; j < size_y ; j++)
  110.         for(i = 0 ; i < (((size_x+2)/2)*2) ; i ++) {
  111.             x = pa[i+j*ld1+k*ld1*ld2] - pb[i+j*ld1+k*ld1*ld2];
  112.             if( (t= x*x) > (emax)) {
  113.             nerr++;
  114.             }
  115.             err += t;
  116.             if( t > maxerr) maxerr = t;
  117.         }
  118.     if( nerr == 0 ){
  119.         fprintf( stderr, ".");
  120.     }
  121.     else { 
  122.         fprintf( stderr, "@%d@", nerr);
  123.     }
  124.     inc = 1;
  125.     x = 0;
  126.     for( j= 0, ptr = pRef ; j < size_z ; j++, ptr += ld1*ld2) { 
  127.         for ( i = 0, pt = ptr; i < size_y; i++, pt += ld1) { 
  128.         x +=  THIS_SUM( &size_x, pt, &inc);
  129.         }
  130.     }
  131.     dx = x - pa[0];
  132.     if( (t= dx *dx) > TOLERANCE){
  133.         fprintf( stderr, 
  134.         "\nError direct FFT(%d,%d) at index #0 : (%f)<->(%f)error(%f)",
  135.         size_x, size_y, pa[0], x, dx);
  136.     }
  137.     dx = x - pb[0];
  138.     if( (t= dx *dx) > TOLERANCE){
  139.         fprintf( stderr, 
  140.         "\nError direct DFT(%d,%d) at index #0 : (%f)<->(%f)error(%f)",
  141.         size_x, size_y, pb[0],x,dx);
  142.     }
  143. /* *******************************************************
  144.     NOT SO PACKED inverse Fourier Transform
  145. ******************************************************* */
  146.     is_wrong = THIS_FFTU( 1, size_x, size_y, size_z, pa, ld1, ld2, pSave);
  147.  
  148.     t = 1. / (double)(size_x * size_y * size_z);
  149.  
  150.     THIS_SCAL( size_x, size_y, size_z, t, pa, ld1, ld2);
  151.  
  152.     emax = TOLERANCE;
  153.     emax = emax * emax;
  154.  
  155.     sum = 0.;
  156.     err= 0.;
  157.     nerr= 0;
  158.     maxerr= 0.;
  159.     for( k = 0 ; k < size_z ; k++)
  160.         for( j = 0 ; j < size_y ; j++)
  161.         for(i = 0 ; i < size_x ; i ++) {
  162.             sum += (pRef[i+j*ld1+k*ld1*ld2] * pRef[i+j*ld1+k*ld1*ld2]);
  163.             x = pa[i+j*ld1+k*ld1*ld2] - pRef[i+j*ld1+k*ld1*ld2];
  164.             if( (t= x*x) > (emax))
  165.             nerr++;
  166.             err += t;
  167.             if( t > maxerr) maxerr = t;
  168.     }
  169.     err = sqrt(err / (double)(size_x*size_y));
  170.     sum = sqrt(sum / (double)(size_x*size_y));
  171.     maxerr = sqrt(maxerr);
  172.     printf(": avg:(%8.3g /%8.3g)= %8.3g <> maxerr: (%8.3g /%8.3g)= %8.3g", 
  173.     err, sum, err/sum, maxerr, sum, maxerr/sum);
  174.     if(nerr){
  175.         printf("\n%d errors detected \n", nerr);
  176.         exit(-2);
  177.     }
  178.  
  179.     free(pa);
  180.     }
  181.     printf("\n");
  182.     return(0);
  183. }
  184.  
  185.  
  186. int is_random;
  187.  
  188. void GetArguments( argc, argv)
  189. int argc;
  190. char *argv[];
  191. {
  192.     int i, j, k;
  193.     int nerror = 0;
  194.  
  195. #define ON    1
  196.  
  197.     n_trials = DEF_TRIALS;
  198.     is_random = 1;
  199.  
  200. /* ******************************************************* */
  201.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  202.     if( argv[i][0] == '-') {
  203.         switch ( argv[i][1]) {
  204.         case 'i' :
  205.         case 'I' :  
  206.                 is_random = 0;
  207.                 break;
  208.         default  :  nerror = ON;
  209.         }
  210.     }
  211.     else {
  212.         if( i+1 > argc)
  213.         nerror = ON;
  214.         else { 
  215.         n_trials = atoi( argv[i]);
  216.         }
  217.     }
  218.     }
  219.     if( nerror == ON) {
  220.     fprintf( stderr, 
  221.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  222.     exit(-1);
  223.     }
  224. }
  225.  
  226. void get_values()
  227. {
  228.   if( is_random){ 
  229.     size_x = random() % MAX_3D + 1;
  230.     size_y = random() % MAX_3D + 1;
  231.     size_z = random() % MAX_3D + 1;
  232.     if( !(n_trials % 5))
  233.     size_z = size_y = size_x;
  234.   } else{
  235.     printf( "Enter SIZE_X : ");
  236.     scanf("%d", &size_x);
  237.     printf( "Enter SIZE_Y : ");
  238.     scanf("%d", &size_y);
  239.     printf( "Enter SIZE_Z : ");
  240.     scanf("%d", &size_z);
  241.   }
  242. }
  243.